##########################################################################################

library(data.table)
library(deconstructSigs)
library(SignatureEstimation)
library(ggplot2)
library(RColorBrewer)
library(reshape2)
library(pheatmap)
library(optparse)
library(ggalluvial)
library(ggsci)

##########################################################################################

option_list <- list(
    make_option(c("--maf_cancer_mss_file"), type = "character") ,
    make_option(c("--maf_cancer_msi_file"), type = "character") ,
    make_option(c("--maf_im_mss_file"), type = "character") ,
    make_option(c("--maf_im_msi_file"), type = "character") ,
    make_option(c("--cosmic_sig_file"), type = "character") ,
    make_option(c("--out_path"), type = "character") ,
    make_option(c("--info_file"), type = "character")
)

if(1!=1){
    
    maf_cancer_mss_file <- "~/20220915_gastric_multiple/dna_combinePublic/maf_public/All_use.maf"
    maf_cancer_msi_file <- "~/20220915_gastric_multiple/dna_combinePublic/maf_public/All_use.msi.maf"
    
    maf_im_mss_file <- "~/20220915_gastric_multiple/dna_combinePublic/maf_public/All_use.IM.maf"
    maf_im_msi_file <- "~/20220915_gastric_multiple/dna_combinePublic/maf_public/All_use.IM.msi.maf"

    info_file <- "~/20220915_gastric_multiple/dna_combinePublic/public_ref/combine/MutationInfo.combine.addMolecularSubType.rmMIX.tsv"

    cosmic_sig_file <- "~/ref/MutationSignature/COSMIC_v3.4_SBS_GRCh37.txt"

    out_path <- "~/20220915_gastric_multiple/dna_combinePublic/finalPlot/revise/signature"

}

###########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

maf_cancer_mss_file <- opt$maf_cancer_mss_file
maf_cancer_msi_file <- opt$maf_cancer_msi_file
maf_im_mss_file <- opt$maf_im_mss_file
maf_im_msi_file <- opt$maf_im_msi_file
cosmic_sig_file <- opt$cosmic_sig_file
info_file <- opt$info_file
out_path <- opt$out_path
dir.create(out_path , recursive = T)

###########################################################################################

info <- data.frame(fread(info_file))
dat_maf_im_mss <- data.frame(fread( maf_im_mss_file ))
dat_maf_im_msi <- data.frame(fread( maf_im_msi_file ))
dat_maf_gc_mss <- data.frame(fread( maf_cancer_mss_file ))
dat_maf_gc_msi <- data.frame(fread( maf_cancer_msi_file ))

###########################################################################################
## 修改突变信号格式
dat_cosmic_sig <- data.frame(fread( cosmic_sig_file ))
rownames(dat_cosmic_sig) <- dat_cosmic_sig$Type
dat_cosmic_sig <- dat_cosmic_sig[,-1]

false_sbs <- c("SBS27" , "SBS43" , "SBS45" , 
    "SBS46" , "SBS47" , "SBS48" , 
    "SBS49" , "SBS50" , "SBS51" , 
    "SBS52" , "SBS53" , "SBS54" , 
    "SBS55" , "SBS56" , "SBS57" , 
    "SBS58" , "SBS59" , "SBS60" , 
    "SBS95"
    )

dat_cosmic_sig <- dat_cosmic_sig[,!colnames(dat_cosmic_sig) %in% false_sbs]

###########################################################################################
## 合并数据
dat_combine_gc <- rbind(dat_maf_gc_mss , dat_maf_im_msi)
dat_combine_im <- rbind(dat_maf_im_mss , dat_maf_im_msi)

###########################################################################################
## 计算信号
computeMutsig <- function( dat_use = dat_use){

    sigs.input <- mut.to.sigs.input(dat_use, sample.id = "type", chr = "Chromosome", pos = "Start_position",
                                        ref = "Reference_Allele", alt = "Tumor_Seq_Allele2")
    sigs.input_plot <- sigs.input
    Sig.input.t <- t(sigs.input)
    sigs.input <- apply(sigs.input, 1, function(x) as.numeric(x))
    sigs.input.N.all <- apply(sigs.input, 1, function(x) x/sum(x))
    sigs.input.N.all <- as.data.frame(t(sigs.input.N.all))

    select <- rownames(dat_cosmic_sig)
    input <- Sig.input.t[match(select,rownames(Sig.input.t)),]
    E1 = findSigExposures(input, dat_cosmic_sig, decomposeQP)
    res <- E1$exposures

    return(res)
}

##############################################
## 分析不同分子亚型的IM，其突变信号特征
info_njmu <- subset(info , From=="NJMU")
dat_use <- merge( dat_combine_im , info_njmu[,c("Tumor" , "Molecular.subtype")] , by.x = "Tumor_Sample_Barcode" , by.y = "Tumor" )
dat_use2 <- subset(dat_use , Molecular.subtype != "MSI")
dat_use2$Molecular.subtype <- "MSS"
dat_use <- rbind(dat_use , dat_use2)
dat_use$type <- paste0("IM_" , dat_use$Molecular.subtype)
mutsig_im_mol <- computeMutsig( dat_use = dat_use)
out_file <- paste0( out_path , "/SBS96_deconstuctSig.IM_moltype.csv" )
write.csv( mutsig_im_mol , out_file )

##############################################
## 在肠化中分饮酒与非饮酒，比较突变信号差异
info_njmu <- subset(info , From == "NJMU" & Alcohol != "unknown" )
dat_use <- merge( dat_combine_im , info_njmu[,c("Tumor", "Molecular.subtype" , "Alcohol")] , by.x = "Tumor_Sample_Barcode" , by.y = "Tumor" )
dat_use$Molecular.subtype <- ifelse( dat_use$Molecular.subtype == "GS" , "GS" , "Other" )
dat_use$type <- paste0("IM_" , dat_use$Molecular.subtype , "_" , dat_use$Alcohol)
mutsig_im_Alcohol <- computeMutsig( dat_use = dat_use)
out_file <- paste0( out_path , "/SBS96_deconstuctSig.IM_drink.csv" )
write.csv( mutsig_im_Alcohol , out_file )

##############################################
## 在肠化中分吸烟与非吸烟，比较突变信号差异
info_njmu <- subset(info , From == "NJMU" & MS_Type!="MSI" & Tobacco != "unknown")
dat_use <- merge( dat_combine_im , info_njmu[,c("Tumor" , "Tobacco")] , by.x = "Tumor_Sample_Barcode" , by.y = "Tumor" )
dat_use$type <- paste0("IM_" , dat_use$Tobacco)
mutsig_im_smoke <- computeMutsig( dat_use = dat_use)
out_file <- paste0( out_path , "/SBS96_deconstuctSig.IM_smoke.csv" )
write.csv( mutsig_im_smoke , out_file )

##############################################
## 比较IGC和DGC的突变信号差异，分MSS、GS、CIN和MSI
info_mss <- subset( info , MS_Type!="MSI" )
info_gs <- subset( info , Molecular.subtype!="GS" )
info_cin <- subset( info , Molecular.subtype!="CIN" )
info_msi <- subset( info , MS_Type=="MSI" )
dat_use_mss <- merge( dat_combine_gc , info_mss[,c("Tumor" , "Class")] , by.x = "Tumor_Sample_Barcode" , by.y = "Tumor" )
dat_use_mss$type <- paste0( dat_use_mss$Class , "_MSS" )

dat_use_gs <- merge( dat_combine_gc , info_gs[,c("Tumor" , "Class")] , by.x = "Tumor_Sample_Barcode" , by.y = "Tumor" )
dat_use_gs$type <- paste0( dat_use_gs$Class , "_GS" )

dat_use_cin <- merge( dat_combine_gc , info_gs[,c("Tumor" , "Class")] , by.x = "Tumor_Sample_Barcode" , by.y = "Tumor" )
dat_use_cin$type <- paste0( dat_use_cin$Class , "_CIN" )

dat_use_msi <- merge( dat_combine_gc , info_msi[,c("Tumor" , "Class")] , by.x = "Tumor_Sample_Barcode" , by.y = "Tumor" )
dat_use_msi$type <- paste0( dat_use_msi$Class , "_MSI" )

dat_use <- rbind(dat_use_mss , dat_use_gs , dat_use_cin , dat_use_msi)
mutsig_gc_mol <- computeMutsig( dat_use = dat_use)

out_file <- paste0( out_path , "/SBS96_deconstuctSig.GC_moltype.csv" )
write.csv( mutsig_gc_mol , out_file )

##############################################
## 比较IGC和DGC的突变信号差异，分饮酒和不饮酒
info_mss_drink <- subset( info , MS_Type!="MSI" & Alcohol != "unknown" )
dat_use <- merge( dat_combine_gc , info_mss_drink[,c("Tumor" , "Class" , "Alcohol")] , by.x = "Tumor_Sample_Barcode" , by.y = "Tumor" )
dat_use$type <- paste0( dat_use$Class , "_" , dat_use$Alcohol)

mutsig_gc_drink <- computeMutsig( dat_use = dat_use)

out_file <- paste0( out_path , "/SBS96_deconstuctSig.MSSGC_drink.csv" )
write.csv( mutsig_gc_drink , out_file )

##############################################
## 比较IGC和DGC的突变信号差异，分吸烟和不吸烟
info_mss_smoke <- subset( info , MS_Type!="MSI" & Tobacco != "unknown" )
dat_use <- merge( dat_combine_gc , info_mss_smoke[,c("Tumor" , "Class" , "Tobacco")] , by.x = "Tumor_Sample_Barcode" , by.y = "Tumor" )
dat_use$type <- paste0( dat_use$Class , "_" , dat_use$Tobacco)

mutsig_gc_smoke <- computeMutsig( dat_use = dat_use)

out_file <- paste0( out_path , "/SBS96_deconstuctSig.MSSGC_smoke.csv" )
write.csv( mutsig_gc_smoke , out_file )


###########################################################################################
## 提取关注的信号占比
sbs_intersting <- c("SBS4" , "SBS16" , "SBS22a" , "SBS24")

plotFun <- function(tmp_dat = tmp_dat , out_file = out_file , width =  width){
    all_class <- colnames(tmp_dat)
    sig_s <- c()
    for(j in 1:length(all_class)){
        sig_s <- c(sig_s,c(names(which(tmp_dat[,all_class[j]]>0.02))))
    }

    ALL_sig <- unique(c(sig_s , sbs_intersting))
    all_sig <- c()
    for(j in 1:length(all_class)){
        all_sig <- rbind(all_sig,data.frame(
            Class = all_class[j] ,
            Exposures = c(tmp_dat[ALL_sig,all_class[j]],
            Others=1-sum(tmp_dat[ALL_sig,all_class[j]])) ,
            Sig = c(ALL_sig,"Others") ))
    }

    #堆叠柱状图
    all_sig$Class_x <- sapply(strsplit(as.character(all_sig$Class) , "_") , "[" , 2)
    all_sig$Class <- factor(all_sig$Class, levels=unique(all_sig$Class)[order(unique(all_sig$Class))], ordered=TRUE) 
    all_sig$Sig <- factor(all_sig$Sig, levels=unique(all_sig$Sig), ordered=TRUE) 
    all_sig$type <- sapply(strsplit(as.character(all_sig$Class) , "_") , "[" , 1)

    p <- ggplot(all_sig, aes(x = Class_x, y= Exposures, fill = Sig,stratum=Sig, alluvium=Sig)) + 
      geom_col(width = 0.5, color='black')+
      facet_grid(.~type , scales = "free" ) +
      geom_flow(width=0.5,alpha=0.4, knot.pos=0)+
      labs(x='Group',y = 'Ratio')+
      scale_fill_manual(values = pal_ucscgb("default", alpha = 0.6)(26)) +
      theme(
            legend.position = 'right',
            legend.title = element_blank() ,
            panel.grid.major=element_blank(),
            panel.grid.minor=element_blank(),
            panel.background = element_blank(),
            panel.border = element_blank(),
            plot.title = element_text(size = 12,color="black",face='bold'),
            legend.text = element_text(size = 12,color="black",face='bold'),
            axis.text.y = element_text(size = 12,color="black",face='bold'),
            axis.title.x = element_text(size = 12,color="black",face='bold'),
            axis.title.y = element_text(size = 12,color="black",face='bold'),
            axis.text.x = element_text(size = 12,color="black",face='bold') ,
            axis.ticks.length = unit(0.2, "cm") ,
            strip.text.x = element_text(size = 15, colour = "black",face='bold') ,
            axis.line = element_line(size = 0.5)) 

    ggsave(out_file, p, width = width , height = 6.5)
}

## 提取感兴趣的信号以及在任一分类中占比超过0.05的信号
tmp_dat <- mutsig_gc_drink
out_file <- paste0( out_path , "/SBS96_deconstuctSig.MSSGC_drink.plot.pdf" )
width <- 8
plotFun(tmp_dat = tmp_dat , out_file = out_file , width =  width)

tmp_dat <- mutsig_gc_smoke
out_file <- paste0( out_path , "/SBS96_deconstuctSig.MSSGC_smoke.plot.pdf" )
width <- 7.5
plotFun(tmp_dat = tmp_dat , out_file = out_file , width =  width)

tmp_dat <- mutsig_gc_smoke
out_file <- paste0( out_path , "/SBS96_deconstuctSig.MSSGC_smoke.plot.pdf" )
width <- 7.5
plotFun(tmp_dat = tmp_dat , out_file = out_file , width =  width)

tmp_dat <- mutsig_im_Alcohol
out_file <- paste0( out_path , "/SBS96_deconstuctSig.IM_drink.plot.pdf" )
width <- 4.5
plotFun(tmp_dat = tmp_dat , out_file = out_file , width =  width)

tmp_dat <- mutsig_im_smoke
out_file <- paste0( out_path , "/SBS96_deconstuctSig.IM_smoke.plot.pdf" )
width <- 4.5
plotFun(tmp_dat = tmp_dat , out_file = out_file , width =  width)